library(foreign)
library(AER)

source("joint-iv-utils.R")


## Code from Hansen & Bowers (2009)
hh <- read.dta("NHrep_household.dta", convert.underscore=TRUE)
hh$AGE1[hh$AGE1MISS==1] <- NA
hh$AGE2[hh$AGE2MISS==1] <- NA
hh$V98.1[hh$V98.1==99] <- NA
hh$V98.2[hh$V98.2==99] <- NA
hh$V96.2.0[hh$PERSONS==1] <- hh$V96.1.0[hh$PERSONS==1]
hh$hh.n <-   hh$PERSONS
if (!('ivdl' %in% ls()))
   ivdl <- read.dta('NHrep_individual.dta',convert.underscore=TRUE)
hh$hh.ty <- ifelse(hh$PERSONS==1, ifelse(is.na(hh$V98.1),0,hh$V98.1),
                   ifelse(is.na(hh$V98.1),0,hh$V98.1) +
                   ifelse(is.na(hh$V98.2),0,hh$V98.2) )
hh$ipc <- as.logical(ivdl$cntany[match(hh$ID1, ivdl$id1)])
hh$phc <- as.logical(ivdl$pcntany[match(hh$ID1, ivdl$id1)])
hh$mailings <- ivdl$mailings[match(hh$ID1, ivdl$id1)]
hh$phongotv <- as.logical(hh$PHONGOTV)
hh$persngrp <- as.logical(hh$PERSNGRP)
hh$mailgrp <- as.logical(hh$MAILGRP)
row.names(hh) <- as.character(hh$ID1)
hh$ccs <- complete.cases(hh$hh.ty)
hh$phonegrp <- hh$PHONGOTV | hh$BLOOD!="not selected"
hh$any.to <- 1 * (hh$hh.ty > 0)
hh.mail <- hh[hh$mailgrp,]
hh.nomail <- hh[!hh$mailgrp,]
hh <- hh.nomail

nrow(hh)

## treatment exclusion test
phn.tet <- lm(phc ~ persngrp, data = hh, subset = phonegrp == 1)
summary(phn.tet)
confint(phn.tet)
ipc.tet <- lm(ipc ~ phonegrp, data = hh, subset = persngrp == 1)
summary(ipc.tet)
confint(ipc.tet)

phn.out <- ivreg(hh.ty ~ phc | phonegrp, data = hh)
summary(phn.out)
confint(phn.out)

person.out <- ivreg(hh.ty ~ ipc | persngrp, data = hh)
summary(person.out)
confint(person.out)

out.itt <- lm(hh.ty ~ phonegrp*persngrp, data = hh)
summary(out.itt)


out <- ivreg(hh.ty ~ phc*ipc | phonegrp*persngrp, data = hh)
summary(out)
confint(out, level = 0.95)
confint(out, level = 0.90)

out <- ivreg(hh.ty ~ phc+ipc  | phonegrp+persngrp, data = hh)
summary(out)
confint(out, level = 0.95)
confint(out, level = 0.90)

out.inter <- ivreg(hh.ty ~ phc*ipc | phonegrp*persngrp, data = hh)
summary(out.inter)
joint.eff <- sum(coef(out.inter)[c(2,3,4)])
tmp <- vcov(out.inter)
this.se <- tmp[2,2] + tmp[3,3] + tmp[4,4] + 2*tmp[2,3] + 2*tmp[2,4] + 2*tmp[3,4]
joint.eff
joint.eff + qt(0.975, df = out.inter$df.residual)*this.se
joint.eff - qt(0.975, df = out.inter$df.residual)*this.se
prsn.eff1 <- sum(coef(out.inter)[c(3,4)])
prsn.eff1.se <- tmp[3,3] + tmp[4,4] + 2*tmp[3,4]
prsn.eff1 + qt(0.975, df = out.inter$df.residual)*prsn.eff1.se
prsn.eff1 - qt(0.975, df = out.inter$df.residual)*prsn.eff1.se
phn.eff1 <- sum(coef(out.inter)[c(2,4)])
phn.eff1.se <- tmp[2,2] + tmp[4,4] + 2*tmp[2,4]
phn.eff1 + qt(0.975, df = out.inter$df.residual)*phn.eff1.se
phn.eff1 - qt(0.975, df = out.inter$df.residual)*phn.eff1.se


out <- joint.iv(y=hh$hh.ty, d1 = hh$ipc, d2 = hh$phc, z1 = hh$persngrp, z2 = hh$phonegrp)
out$rho
laje.est <- out$tau["laje"]
laje.ci95 <- laje.est + qnorm(c(0.025, .975)) * out$se["laje"]
laje.ci90 <- laje.est + qnorm(c(0.05, .95)) * out$se["laje"]

laie.est <- out$tau["laie"]
laie.ci95 <- laie.est + qnorm(c(0.025, .975)) * out$se["laie"]
laie.ci90 <- laie.est + qnorm(c(0.05, .95)) * out$se["laie"]

cc.10.est <- out$tau["cc.10"]
cc.10.ci95 <- cc.10.est + qnorm(c(0.025, .975)) * out$se["cc.10"]
cc.10.ci90 <- cc.10.est + qnorm(c(0.05, .95)) * out$se["cc.10"]

cc.11.est <- out$tau["cc.11"]
cc.11.ci95 <- cc.11.est + qnorm(c(0.025, .975)) * out$se["cc.11"]
cc.11.ci90 <- cc.11.est + qnorm(c(0.05, .95)) * out$se["cc.11"]

cc.20.est <- out$tau["cc.20"]
cc.20.ci95 <- cc.20.est + qnorm(c(0.025, .975)) * out$se["cc.20"]
cc.20.ci90 <- cc.20.est + qnorm(c(0.05, .95)) * out$se["cc.20"]

cc.21.est <- out$tau["cc.21"]
cc.21.ci95 <- cc.21.est + qnorm(c(0.025, .975)) * out$se["cc.21"]
cc.21.ci90 <- cc.21.est + qnorm(c(0.05, .95)) * out$se["cc.21"]

cn.est <- out$tau["cn"]
cn.ci95 <- cn.est + qnorm(c(0.025, .975)) * out$se["cn"]
cn.ci90 <- cn.est + qnorm(c(0.05, .95)) * out$se["cn"]

nc.est <- out$tau["nc"]
nc.ci95 <- nc.est + qnorm(c(0.025, .975)) * out$se["nc"]
nc.ci90 <- nc.est + qnorm(c(0.05, .95)) * out$se["nc"]

cc.cn.est <- out$tau["cc.cn"]
cc.cn.ci95 <- cc.cn.est + qnorm(c(0.025, .975)) * out$se["cc.cn"]
cc.cn.ci90 <- cc.cn.est + qnorm(c(0.05, .95)) * out$se["cc.cn"]

cc.nc.est <- out$tau["cc.nc"]
cc.nc.ci95 <- cc.nc.est + qnorm(c(0.025, .975)) * out$se["cc.nc"]
cc.nc.ci90 <- cc.nc.est + qnorm(c(0.05, .95)) * out$se["cc.nc"]


ipc.pts <- c(coef(person.out)[2], coef(out.inter)[3], prsn.eff1, cc.10.est, cn.est, cc.cn.est, cc.11.est)
ipc.ci95.hi <- c(confint(person.out)[2,2],confint(out.inter)[3,2], prsn.eff1 + qt(0.975, df = out.inter$df.residual)*prsn.eff1.se, cc.10.ci95[2], cn.ci95[2], cc.cn.ci95[2], cc.11.ci95[2])
ipc.ci95.lo <- c(confint(person.out)[2,1],confint(out.inter)[3,1], prsn.eff1 - qt(0.975, df = out.inter$df.residual)*prsn.eff1.se, cc.10.ci95[1], cn.ci95[1], cc.cn.ci95[1], cc.11.ci95[1])
ipc.ci90.hi <- c(confint(person.out, level = 0.9)[2,2],confint(out.inter, level = 0.9)[3,2], prsn.eff1 + qt(0.95, df = out.inter$df.residual)*prsn.eff1.se, cc.10.ci90[2], cn.ci90[2], cc.cn.ci90[2], cc.11.ci90[2])
ipc.ci90.lo <- c(confint(person.out, level = 0.9)[2,1],confint(out.inter, level = 0.9)[3,1], prsn.eff1 - qt(0.95, df = out.inter$df.residual)*prsn.eff1.se, cc.10.ci90[1], cn.ci90[1], cc.cn.ci90[1], cc.11.ci90[1])
ipc.labs <- c("LATE","iTSLS (Phone = 0)", "iTSLS (Phone = 1)", expression(tau["cc,2"](0)), expression(tau["nc,2"](0)), expression(tau["-c,2"](0)), expression(tau["cc,2"](1)))


phn.pts <- c(coef(phn.out)[2], coef(out.inter)[2], phn.eff1, cc.20.est, nc.est, cc.nc.est, cc.21.est)
phn.ci95.hi <- c(confint(phn.out)[2,2],confint(out.inter)[2,2], phn.eff1 + qt(0.975, df = out.inter$df.residual)*phn.eff1.se, cc.20.ci95[2], nc.ci95[2], cc.nc.ci95[2], cc.21.ci95[2])
phn.ci95.lo <- c(confint(phn.out)[2,1],confint(out.inter)[2,1], phn.eff1 - qt(0.975, df = out.inter$df.residual)*phn.eff1.se, cc.20.ci95[1], nc.ci95[1], cc.nc.ci95[1], cc.21.ci95[1])
phn.ci90.hi <- c(confint(phn.out, level = 0.9)[2,2],confint(out.inter, level = 0.9)[2,2], phn.eff1 + qt(0.95, df = out.inter$df.residual)*phn.eff1.se, cc.20.ci90[2], nc.ci90[2], cc.nc.ci90[2], cc.21.ci90[2])
phn.ci90.lo <- c(confint(phn.out, level = 0.9)[2,1],confint(out.inter, level = 0.9)[2,1], phn.eff1 - qt(0.95, df = out.inter$df.residual)*phn.eff1.se, cc.20.ci90[1], nc.ci90[1], cc.nc.ci90[1], cc.21.ci90[1])

pts <- c(joint.eff, coef(out.inter)[4], laje.est, laie.est)
ci95.hi <- c(joint.eff + qt(0.975, df = out.inter$df.residual)*this.se, confint(out.inter)[4,2], laje.ci95[2], laie.ci95[2])
ci95.lo <- c(joint.eff - qt(0.975, df = out.inter$df.residual)*this.se, confint(out.inter)[4,1], laje.ci95[1], laie.ci95[1])
ci90.hi <- c(joint.eff + qt(0.95, df = out.inter$df.residual)*this.se, confint(out.inter, level = 0.9)[4,2], laje.ci90[2], laie.ci90[2])
ci90.lo <- c(joint.eff - qt(0.95, df = out.inter$df.residual)*this.se, confint(out.inter, level = 0.9)[4,1],  laje.ci90[1], laie.ci90[1])


y.pts <- c(1,4,2,5)
col.pts <- c("indianred", "indianred", "dodgerblue", "dodgerblue")

y.pts2 <- c(-1,1,6,3,4,2,7)
these.cols <- c("dodgerblue", "indianred", "indianred", "dodgerblue", "dodgerblue", "dodgerblue", "dodgerblue")
#cairo_pdf(file = "../figs/gerber-green.pdf",width = 10, height = 4, pointsize = 16, family = "Minion Pro")
par(mfrow = c(1,3), mar = c(4.1, 0.1, 2.1, 0.1))
plot(x = NA, y = NA, xlim = c(-3,3), ylim = c(-1.5,7.5), las = 1, xlab = "Treatment Effects", ylab="", bty = "n", yaxt = "n", main = "Phone Contact", xaxt = "n")
axis(side = 1, at = c(-2,-1, 0, 1, 2))
points(x = phn.pts, y = y.pts2, pch = 19, col = these.cols)
segments(x0 = phn.ci95.lo, x1 = phn.ci95.hi, y0 = y.pts2, col = these.cols)
segments(x0 = phn.ci90.lo, x1 = phn.ci90.hi, y0 = y.pts2, lwd = 3, col = these.cols)
abline(v=0, col = "grey70", lty = 2)
text(x = -3.25, y = y.pts2, c("LATE","iTSLS (In-Person = 0)", "iTSLS (In-Person = 1)", expression(tau["cc,1"](0)), expression(tau["cn,1"](0)), expression(tau["c-,1"](0)), expression(tau["cc,1"](1))), pos = 4, col = these.cols, cex = c(1, 1, 1, rep(1.5, 4)))
plot(x = NA, y = NA, xlim = c(-3,3), ylim = c(-1.5,7.5), las = 1, xlab = "Treatment Effects", ylab="", bty = "n", yaxt = "n", main = "In-Person Contact", xaxt = "n")
axis(side = 1, at = c(-2,-1, 0, 1, 2))
points(x = ipc.pts, y = y.pts2, col = these.cols, pch = 19)
segments(x0 = ipc.ci95.lo, x1 = ipc.ci95.hi, y0 = y.pts2, col = these.cols)
segments(x0 = ipc.ci90.lo, x1 = ipc.ci90.hi, y0 = y.pts2, lwd = 3, col = these.cols)
abline(v=0, col = "grey70", lty = 2)
text(x = -3.25, y = y.pts2, ipc.labs, pos = 4, col = these.cols, cex = c(1, 1, 1, rep(1.5, 4)))
plot(x = NA, y = NA, xlim = c(-3,3), ylim = c(-1.5,5.5), las = 1, xlab = "Treatment Effects", ylab="", bty = "n", yaxt = "n", main = "Joint and Interaction Effects", xaxt = "n")
axis(side = 1, at = c(-2,-1, 0, 1, 2))
points(x = pts, y = y.pts, pch = 19, col = col.pts)
segments(x0 = ci95.lo, x1 = ci95.hi, y0 = y.pts, col = col.pts)
segments(x0 = ci90.lo, x1 = ci90.hi, y0 = y.pts, lwd = 3, col = col.pts)
abline(v=0, col = "grey70", lty = 2)
text(x = -3.25, y = y.pts, c("iTSLS Joint", "iTSLS Interaction", expression(tau["LAJE"]), expression(tau["LAIE"])), pos = 4, col = col.pts, cex = c(1,1, rep(1.5, 2)))
#dev.off()
